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ABSTRACT 



In this paper I review and discuss the basic concepts of accretion disks, focused especially on 

fr^ , the case of accretion disks around black holes. The well known a-model is revisited, showing 

the strengths and weaknesses of the model. Other turbulent viscosity prescription, based 

ftV' on the Reynolds number, that may improve our understanding of the accretion paradigm is 

hH . discussed. A simple but efficient mathematical model of a self-gravitating accretion disk, as 

• I well as observational evidence of these objects, are also included. 

O ' Subject headings: supermassive black holes - accretion disks - active galactic nuclei 



00 

rn 

O 



1. Brief review of accretion disk theory 

Accretion disks are expected to be associated to dif- 
ferent astrophysical objects in the universe and, among 
them it is possible to mention; disks around young stars 
(T Tauri), probably sites of planetary formation; disks in 
close binary systems like /3-Lyra, cataclysmic variables, 
novae. X-ray binaries, etc., and AGNs in general. 

The basic idea of accretion disks was proposed soon 
after the confirmation of quasars, and its classical picture 
established at the beginning of the 1970s. 

When the gas infalls into the compact object, some of 
the binding gravitational energy of the gas is released due 
the viscosity present in the accretion flow. The amount 
of gravitational energy liberated during the process is 
huge. This was first pointed out by Zeldovich (1964), Zel- 
dovich & Novikov (1967), and independently by Salpeter 
(1964). However, they considered an isolated compact 
object converting gravitational energy into radiation dur- 
ing its collapse. 

Consequently to these historical papers, several mod- 
els explaining observed energetic phenomena like X-ray 
binaries appeared. For instance Shklovsky (1967) ex- 
plained X-ray emissions from Sco X-1^ as accretion of 
matter from a companion onto a neutron star. However 
at the time, the majority of investigators assumed in- 
stead Bremsstrahlung radiation of an optically thin layer 
of hot plasma to explain energetic objects. 

Under the same argumentation that gravitational en- 
ergy can be converted into radiation by viscous forces, 
Lyndcn-Bell proposed in 1969 that the energy source of 
quasars can be explained by viscous accretion disks. 

The accretion disk model grew in popularity and be- 
gan to be recognized as the theory capable to describe 
important energetic emissions observed in the Universe. 
It was in 1973 that Shakura & Sunyaev proposed a funda- 
mental theory of accretion disks known as the standard 
accretion disk model or simply the a-disk model. Very 
soon, in the same year, it was generalized to the rela- 
tivistic version by Novikov & Thornc (1973). 

The Shakura & Sunyaev model, described a relatively 
cold disk (in range ^ 10^ — 10'"' K for typical case of 
quasars) due an efficient process of conversion and evac- 
uation of gravitational energy into radiation. The disk 
was also assumed geometrically thin and optically thick. 
They introduced a clever parameterization for the turbu- 
lent viscosity, needed to guarantee the transfer of angular 
momentum and the process of accretion flow. 

Around the 80s different features of accretion disk 



were explored such as hot disks, radiatively inefficient 
disk, geometrically thick disks, etc. 

It is the case, for instance, of the Polish doughnut 
disk, consisting in an optically thick torus introduced by 
Paczyhski & Wiita in 1980, and later, when advection 
were recognized as an important ingredient, the slim disk 
appeared (Abramowicz 1988), giving a new branch of ac- 
cretion disk, where the advective energy transport was in- 
cluded in the energy balance, producing super-Eddington 
accretion rates. 

The importance of these advective models, called 
ADAFs (Advection-Dominated- Accretion-Flows) , grew 
up in the 1990s, and were developed by several au- 
thors (e.g., Narayan & Yi (1994), Abramowicz & Lasota 
(1995)). 

All these models showed a good agreement with ob- 
servations. Accretion disk is now the accepted paradigm 
in order to describe several energetic phenomena either 
on small scales, like disks around Young stellar object 
(YSO), or Cataclysmic variable stars (CV) or on large 
scales, like AGNs or quasars. Table (1), summarizes some 
characteristic properties, such as the mass, size, temper- 
ature, and luminosity of the accretion disk supposed to 
be at the origin of the indicated astrophysical object. 

Table 1: Properties of the accretion disk in various as- 
trophysical objects. 



Central object 


Central object mass 


Disk size 


Disk mass 


cv/ssxs 


~ Mq 


~Re 


< Mq 


XB/BHB 


>3Me 


'-^ il0 


<.Mq 


YSO 


~ Mq 


~ 100 AU 


^Mq 


AGN 


- lO-'' - IO^Mq 


~pc 


~ 10'' - 10^ Mq 



^Sco X-1 was the first extrasolar X-ray source, discovered in 1962 
(Giacconi et al. ) 



The used nomenclature reads- CV: cataclysmic variable, 

SSXS: supersoft X-ray source, XB: X-ray binary, BHB: 

black hole binary, YSO: young stellar object, AGN: active 

galactic nuclei. 

Summarizing, the central idea of accretion disks lies in 
the fact that the accreting masses possess a considerable 
amount of angular momentum per mass unit which has 
to be removed in order to be accreted into the central 
object. What causes the lost of angular momentum is the 
friction caused by turbulent viscosity working between 
adjacent gas layers in the disk. The faster inner layer 
loses angular momentum and infalls slightly, while the 
next (slower) outer layer gains angular momentum, which 
is given away to the next outer layer, and so on, resulting 
in a continuous flow toward the centre whereas angular 
momentum is transported to the outer region (see Figure 
1 for a pictorial view of the accretion process). If the 
disk is massive enough in the external region, the gain 
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Fig. 1. — Schematic diagram of an accretion disk pro- 
cess. The inner edge of the annulus at (r — Sr) receives 
a force —trip{r — Sr) per unit surface, gaining angular 
momentum. Likewise, the outer edge of the annulus at 
r receives a force tr^ per unit surface, giving away an- 
gular momentum. Consequently, angular momentum is 
transfered outward while mass is accreted inward. 



of angular momentum (of the very last outer layer) will 
cause an expansion of the disk because there is no other 
shell of gas to give up the gained extra momenta. In that 
case, material in that layer is transferred outward. 

The friction of the disk heats up the gas, resulting in 
a continuous radiation emission, which is believed to be 
the source of AGNs, quasars. X-ray binaries, luminosities 
(among other phenomena) as mentioned above. 

2. Basic equations 

In order to obtain a realistic model of an accretion 
disk we need to have an adequate description of the flow; 
that is why the complete set of hydrodynamic equations 
must be solved. 

Fluid motion can be completely determined by con- 
servation laws. They state that at all times, a certain 



number of properties, e.g., mass, generalized momentum 
and energy, are conserved quantities. 

We start the analysis with the Navier-Stokes equa- 
tions, describing the motion of a viscous compressible 
fluid with variable dynamic viscosity 77 and assuming bulk 
viscosity coefhcient equal to zero i.e., 773 = 0. 

We consider a fixed volume V wrapped within a sur- 
face 5". The rate of the fluid mass change contained in 
the volume V is 



I / P'^^ 



(1) 



where p is the density of the fluid. 

In the absence of source/sink for matter, Ec^. (1) must 
be equal to the total mass inflow integrated over the sur- 
face. 

Assuming that v is the velocity vector of the fluid, 
the outward mass flow across an element dS is pv ■ dS. 
Therefore, the gained mass by the volume V is calculated 
by integrating over the surface S. Using the divergence 
theorem to express surface integral into volume integral, 
we have, 

- [ pv-dS = - / V • {pv)dV, (2) 

Js Jv 

the minus sign is added because we are considering an 
outflow motion. 

Therefore, equating the above expression to the Eq. 
(1), we flnd. 



d_ 

dt 



[ pdV = - [ V- {pv)dV. 
Jv Jv 



(3) 



From this, the mass conservation law can be written 



as: 



-p + V-(pv) = 0. 



(4) 



This is the usual Eulerian form of the mass continuity 
equation. 

In order to formulate the conservation law for momen- 
tum, it is necessary to define the sources infiuencing the 
variation of momentum. These sources are forces acting 
on the fiuid, they consist in external volume forces fe and 
internal ones fi. The latter depends on the nature of the 
considered fluid, i.e., on the assumptions made about the 
deformations within the fluid and their relation to inter- 
nal stresses. For the later, a viscous stress tensor must 
be introduced, 



nj = ncTij = ^ f [diVj + djVi) - -(V • v)6ij j , (5) 

where aij is the shear tensor and /i the dynamic viscosity 
of the fluid. 

Now, the total internal stresses can be described using 
the stress tensor T. Denoting by P the pressure, we have, 

T = -PI + T, (6) 

where / is the identity matrix, and t the viscous stress 
tensor given by (5). 

Using Equation (6) we can find the internal forces act- 
ing on the fluid, i.e., f ^ = V • T = -VP + V • r. 

If the external forces come from a gravitational poten- 
tial ^, then we have fg = — V^. 

We can now, write down the equation for the momen- 
tum conservation, which reads. 



p(^v+(vV)v 



= VT + p/e 

= -VP + V • T - pV*. (7) 



This equation is commonly referred to as the Navier- 
Stokes equation of motion. For an ideal fluid without 
internal stresses, that is an inviscid flow, i.e., T = 0, this 
reduces to the Euler equation of motion. 

These set of equations are a general description of fluid 
motions. Usually they are very difficult to solve and can- 
not be done without several assumptions about the na- 
ture of the fluid and the geometry of the problem. 

We use cylindrical coordinates (r, (/3, z) to describe our 
disks. From an analytical (and numerical) point of view 
this is a very natural choice, not just due to the disk ge- 
ometry but also because integrations arc performed from 
the last stable orbit riso, and not from the origin of the 
coordinate at r = 0. In that point the Jacobian diverges 
causing several issues in numerical calculations. 

We also take symmetry considerations and vertical in- 
tegration of variables, to reduce the dimension of the 
problem and the number of independent variables. For 
instance, it is assumed an axisymmetric disk, which 
means that all quantities are independent of the (p co- 
ordinate, then ^ = 0. When dropping the azimuthal 
coordinate, it is no longer possible to resolve local dis- 
turbances in the azimuthal direction. This kind of term 
could be a source of viscosity and/or turbulence interest- 
ing to observe, but when we deal with a very extended 
disk (^ 150 pc) computational cost in adding the az- 
imuthal component would be prohibitive. 



As for the vertical integration, the idea is to get rid of 
all z dependencies by integrating the equations through 
the depth of the disk, so we assumed that the flow is 
symmetric with respect to the equatorial plane (mirror 
symmetry about this plane) . This procedure allows us to 
decouple the vertical and radial directions. In that case 
we have for the velocity components Vr, w^, and Vz, 



Vr {r,z)c:i Vr (r) , w^p (r, 2;) ~ v^ (r) , 



0. 



(8) 



The condition Vz — is, of course, equivalent to say 
that the disk is in hydrostatic equilibrium along the z- 
axis. 

Vertical integration is then understood as integration 
over vertical coordinate with the scale of height of the 
disk as limits, i.e., ±iJ(r)/2 . This implies that rather 
than dealing with quantities per unit volume, we will deal 
instead with quantities per unit surface. Integrating the 
density p along the z-axis we obtain the surface density 
given by. 



S(t, r) ^ I pdz 



+H(r)/2 



H(r)/2 



pdz — H{r)p{t,r,z = 0), 

(9) 

where H{r) is the total scale of height. Like the density p, 
all other quantities are vertically averaged. For instance, 
for the viscous stress we obtain, 



T = 



-H{r)/2 



H{r)/2 



"^flV^^' 



(10) 



Using the above symmetry considerations, we can eas- 
ily derive the mass conservation equation, the radial and 
azimuthal components of the Navier-Stokes equation (7) 
within cylindrical coordinates. Therefore, we have: 

2.0.1. mass conservation 

From Eq. (4), we simply obtain, 

as 1 d{rY.Vr 



dt r dr 
2.0.2. the radial component 



0. 



(11) 



The equation for the radial transport, including gas 
pressure, and the three components of viscous stress ten- 
sor T^, (Eqs. (17), (18), (19)), is given by. 



dVr 



dVr 

dr 



1 dP vi 



9* 



S(i,r) dr r dr 



(12) 



Approximations that can be used without strikingly 
changing the physics of the system is to neglect gas gra- 
dient pressure along the radial direction, and assume the 
components T^r — Tr^ = 0, hence oj!'"'^ = 0. These 
assumptions are well established and justified due the 
fact that the azimuthal velocity is greatly superior to the 
radial velocity (i.e., v^ ^ Vr). In that case, the only 
important viscous force lies between two consecutive an- 
nulus r — Ar and r which creates a force per unit surface 
Trip and thus exerts a torque which is responsible for 
transporting the angular momentum. Therefore the last 
Equation (12) simplifies to. 



dvr 
~dt 



dVr 



r 



d^T 



(13) 



dr r 9r ' 

where viscous forces play no role in the radial component. 
The gravitational potential ^t includes the contribution 
from both, the central black hole and the disk itself. Note 
that there is a continued flow of matter as time passes, 
therefore the gravitational potential is also a function 
of time. We then have *r = ^T{r,t) = *b/i(r,i) + 
^disk{r,t). 

2.0.3. the azimuthal component 

The equation for the azimuthal motion includes the 
viscous forces responsible for the angular momentum 
transport. From the equation (7) we have. 






dr 



(14) 



where o^*^"^ and a^'^'^ are the radial and azimuthal accel- 
erations due to the disk viscosity (r and </j-component of 
V • r in cylindrical coordinates), they are given by 



1 



d{rTr, 



rT,{t,r) \ dr 

1 ( djrTr^) 
rJ:(t,r) V dr 



T,„ 



(15) 
(16) 



where the terms Tap of the components of the vertically 
integrated viscous stress tensor (5) arc given in cylindri- 
cal coordinates, by 



the divergence of the velocity field in the above equations, 
can be written as, 



1 5 / N 

V • v= -—-{rvr). 
r or 



(20) 



When V • V = 0, the fiuid is said to be incompressible, 
which is not our case. We let the fluid be compressed 
near the location where forces are applied, that is, its 
density increases locally in response to that force. The 
compressed fluid expands against neighboring fiuid par- 
ticles causing the neighboring fluid itself to compress and 
setting in motion a wave pulse that travels throughout 
the system. 

The radial acceleration due to the viscosity (Eq. 15) 
and the pressure gradient term, present in Equation (12) 
can be neglected, their contributions are negligible next 
to the gravitational forces in the radial direction. Any- 
way, we do take into account these quantities in our cal- 
culations. 

2.0.4. ^^6 scale of height 

The vertical component of Equation (7) can be greatly 
simplified if we assume very small velocities (static sit- 
uation ^ = 0) and no viscous forces in the z-direction 
(i.e., Vz — and Tzz — 0), the disk is being confined to 
the equatorial plane. Then we can neglect all the left- 
hand side term of the equation. The remaining terms 
to balance are the pressure force in the vertical direc- 
tion with the gravitational forces (from both the central 
object and the disk itself). This balance is called the 
hydrostatic equilibrium, and reads: 



IdP 

p dz 



9* 

dz ' 



(21) 



where P should be the total (vertically-averaged) pres- 
sure composed by 



P = P. 



gas 



pCg, the turbulent pressure Ptur = P < 
aT^/3, and a 



where Pgas 

Vf >, the radiation pressure Pad 

7.564 X 10^^'"' erg cm~^ K~^ is the radiation constant. 

The above Eq. (22) can be rewritten as 



T 



T. 



ipip 



Tr, 



2uY.{t,r) 
2vT.{t,r) 
vY,{t, r 



dv. 
dr 3 

Vr 1 

T~ 3 

d /v. 
dr 



V.v 

V-v 



"V 



(17) 
(18) 
(19) 



P 



pcl + p< vf > + 



= PCs 



= PCs 



aT'^ 



< vf > 



aT'^ 
2,pcl 



< vf > 2HaT'^ 



SSc? 



= pclil + e'+^t^H], 



(23) 



where e^ = < w^ >/Cs, is the ratio between the tur- 
bulent velocity and sound speed, and 7 is defined as 
72 = 2aTV3Sc2. 

To calculate the turbulent velocity < Vt > we assume 
an isotropic turbulence. In that case, the kinematical 
turbulent viscosity v is given by 



V = -<Vt><lt>, 



(24) 



where < It > corresponds to the characteristic scale of 
eddies. The turbulent velocity can be approximated us- 
ing < Vt >~< It > ^. Then, the last equation can be 
recast as 



l<Vt>^ 

V = —- 



3 ^ 



(25) 



Matching this expression with the used viscosity pre- 
scription V, and isolating < wt >, we have (using 51 = 



< Vt >2~ 3z^r2. 



(26) 



Therefore, from the definition of the parameter e we 
obtain, 



3un 



(27) 



Defining the total scale of height by the relation 

P 



H 



dp 

dz 



(28) 



and using Equation (21) we can obtain the values oi H{r) 
as a function of the radius r. To find an adequate expres- 
sion, notice that we can always write 



(29) 



dP _ dP dp 
dz dp dz ' 

therefore, using the definition of the effective scale of 
height, we obtain 



H 



dP 1 
dp gz 



(30) 



The first factor in the rhs of the above equation, can 
be obtained from Eq. (23), giving 



^ = cUl + e'+^'H) 



(31) 



and the second one (i.e. l/^z), were gz is the vertical 
component of the gravity acceleration, is taken from g^ — 
-d-^T/dz = GMbhz/r^ + 2i:GY.z/H, in which the total 
gravitational potential was considered. 

With these quantities, Eq. (30) transforms into. 



H 



Pclil 



I'H) 



2itGY. 



nj^H 



(32) 



Resolving the above equation for H, after some alge- 
bra we find a consistent expression for the scale of height 
involving different pressure sources: 



H 



Qn 



K 



(l-/3)-V(l-/3)2+Q (l + e^) 



(33) 

In this equation, there are two important quantities: 
Q and (3. The first one was defined as: 



7=: ^KCs 



ttGS 



(34) 



which gives the limit between a self-gravitating and a ke- 
plerian regime for the scale of height. Given the similar- 
ity, it must not be confounded with the Toomre param- 
eter Q. The latter considers the epicyclic frequency fc, 
while Q considers the Keplerian angular velocity. Only 
in the case of Keplerian motion, when k — U,k, these 
values coincide i.e., Q — Q. 

The second quantity was defined as: 



13: 



aT* 



37rGE2 



(35) 



The parameter /3 measures the relative influence of 
the radiative pressure with respect to the self-gravity of 
the disk. When radiation pressure is increased, it can 
overcome gravity. The disk is then inflated (or it could 
be destroyed), producing a slim disk. 

In this model, we do not impose a thin disk condition 
as in the the case of the Shakura-Sunyaev models where 
H/r << 1. The vertical extent of the disk H is allowed 
to be H/r <; 1. 

Only in the inner region of the disk, where tempera- 
tures are high enough and /3 ^ 1, the scale of height of 
the disk H could be pumped up by radiation pressure 
and be comparable to r, i.e., H ^ r. This force creates 
a small hot corona in this region. Moving away from the 
center, temperatures drop quickly and gravitational tidal 
forces overcome largely the radiation pressure. 

It is important to mention that radiation pressure 
must be taken into account only when the disk is optically 



thick. In the opposite case, radiation passes through the 
gas without exerting any force, hence it is negligible. In 
such a case /? = 0. This condition must to be applied in 
order to calculate the scale of height H from Eq. (33). 

Two interesting limiting cases can be explored from 
the equation for H (Eq. 33). One corresponds to the 
case when the radiation force is negligible (/? ^ 1) and 
self-gravity dominates the entire disk (Q <c 1), applying 
these in Eq. (33), we obtain 



H 



cUl 



27rGS 



(36) 



When the turbulent velocity < w^ > equals the sound 
speed, i.e., e ~ 1, we have a typical expression for the 
scale of height for self-gravitating accretion disk, that 
can be obtained almost directly from the vertical balance 
(Eq. 21), 



H 



ttGS' 



(37) 



The other case arises when the radiation force is still 
unimportant but also, the self-gravity of the disk is neg- 
ligible {Q >• 1), in this case the disk is Keplerian and we 
have. 



H 



,VT 



n 



K 



(38) 



Once the temperature of the disk is calculated from 
the energy equation, it is possible to immediately eval- 
uate the expression for H; recall that scale of height 
is a time-dependent function of the temperature H = 
H{t,T). 

The system of Navicr-Stokes equations still has to 
be supplemented with a specification of the kinematic 
viscosity u (related to the dynamical viscosity through 
v = j-i/ p) as a function of other flow variables. Using 
the same symmetry arguments, vertical integration of the 
kinematic viscosity is defined as. 



The next subsection, the a— disk model is presented 
with their main assumptions, features, successes and lim- 
itations. 

2.1. Standard model: a-disk 

The standard model of accretion disk was first for- 
mulated in the papers of Shakura (1972) and Shakura 
& Sunyaev (1973), and then generalized to the Kerr- 
metric by Novikov & Thorne (1973). It is interesting 
to recall that all these soviet authors (Shakura, Sunyaev, 
Novikov), pioneers in the development of the accretion 
theory, were collaborators of Y.B. Zel'dovich whose in- 
fluence was of extremely importance even if his name 
does not appear in these classical papers. 

The main idea of the a— model was to describe a geo- 
metrically thin non-self gravitating disk by hydrodynam- 
ical equations averaged over the disk thickness. Geomet- 
rically thin and non-self-gravitating disk means that the 
scale of height of the disk H (thickness) is much smaller 
than the radial distance r (i.e., H/r <C 1) and the mass of 
the disk Md is much smaller than the mass of the central 
object M, (i.e., Md ^ M,) so the gravitational influence 
of the disk is negligible. 

The most important process governing the accretion 
of rotating matter is the action of viscous stress within 
the flow. The goal is to describe the motion of a vis- 
cous compressible fluid with variable dynamic viscosity 
rj. Viscous stresses drive accretion by transporting an- 
gular momentum outward and mass inward; it is also a 
way to convert gravitational energy of matter into heat, 
which will be released toward the top and bottom of the 
surface of the disk and radiated away. 

The evolution of disks under viscosity was also ex- 
tensively studied in a paper from Lyndcn-BcU & Pringle 
(1974). 

The original standard disk model (the a-disk model) 
also supposes an optically thick accretion disk and a tur- 
bulent fluid described by a viscous stress tensor which is 
proportional to the total pressure. It is then parametrized 



<v >-- 



1 



-H{r)/2 



p v dz 



(39) 



-H{r)/2 



For simplicity, when we talk about the kinematic vis- 
cosity we mean the z- averaged kinematic viscosity given 
by expression (39). 

The choice of < i^ > is the most obscure and specu- 
lative aspect of accretion theory and a key problem to 
be solved. The main issue is to find a source of r stress 
that could be effective to transport angular momentum 
making accretion possible in a reasonable timescale. 



Trip — apc^ = —aP, (40) 

where a is a dimcnsionlcss constant, that can be fixed 
for values between zero (case when accretion is halted) 
and close to one, and Cg the sound speed of the flow. The 
idea is that viscosity is generated by turbulent motion, 
and the typical turbulence scale It has the same order of 
the vertical disk scale H. On the other hand, turbulent 
motion cannot be supersonic, otherwise energy would be 
dissipated very fast and the typical turbulent velocity 
would drop below the sound speed. 



Using the definition of the viscous tensor Tr^p (Eq. 19), 
the last parameterization (Eq. 40) is completely equiva- 
lent to considerer the kinematical viscosity < v >, given 

by 



< V >— ac^H. 



(41) 



This form of the a— model immediately reflects the 
fact that turbulent viscosity is considered proportional to 
the turbulent velocity vt and the scales of eddies /( in the 
turbulent pattern, i.e., v ex vth, where these quantities 
are supposed to be Vt ~ Cg and It ~ iJ as expressed 
above. Since the turbulent velocity must be Vt <^ Cg and 
the eddies scales It <, H, we obtain two upper limits 
implying that a should be ^ 1. 

It must be clear that the a— model is not a real the- 
ory of viscosity in accretion disk it is just a way to hide 
our ignorance about viscosity, replacing the unknown pa- 
rameter 1/ by another unknown parameter a which is sup- 
posed to be <^ 1. Nonetheless, the model can be easily 
compared with observations constraining even more the 
values of a. For instance, typically quoted values for a 
range from ^ 0.01 for protostellar disks (Hartmann et al. 
1998) to ^ 0.1 for galactic binaries (Lasota 2001). 

The model assumes that the disk is in local thermal 
equilibrium, does not advect heat inwards, and can radi- 
ate the viscous heat efficiently as a blackbody. 

The main assumptions in the standard disk model can 
be summarized as follows: 

• Gravitation is only determined by the central ob- 
ject (no self-gravity effect) 

• The disk is geometrically thin, i.e., H/r <C 1 

• The disk is steady, i.e., d/dt = 

• The disk is axisymmetric i.e., d/dip = 

• Azimuthal motion dominates over radial motion 
i.e., v^ > Vr 

• Hydrostatic balance is assumed in the vertical di- 
rection 

• The disk is optically thick in the vertical direction 

• Viscous heat dissipation balance radiation output 



• The r — If component of the viscous stress tensor 
is assumed proportional to the pressure i.e., Tr^ = 
—aP. All other components are zero 

• There are no magnetic fields present in the disk 



Observations show the electromagnetic energy distri- 
bution in frequency band. The success of an accretion 
disk model lies in the ability to describe, among other 
phenomena, these observations. Of course they depend 
on the nature of the source, whether the object is a bi- 
nary system or an active galaxy for instance. 

This model was quickly recognized; it supplies a ro- 
bust theory of the accretion flow, easy to apply to obser- 
vations. 

It successfully fits, for instance, the emission proper- 
ties of dwarf- novae. X-ray binaries (Cannizzo 1993), UV- 
soft X-ray emission of AGNs (Sun & Malkan 1986). How- 
ever, this picture is unstable under thermal and viscous 
instabilities (this point will be clarified later), moreover, 
other viscosity prescription could work as well in order 
to describe the aforementioned observations. 

2.2. The classical Eddington limit 

The mass M of the central object plays an important 
role in the nature of the accretion disk, especially when 
considering the maximum luminosity allowed by the sys- 
tem. This maximum luminosity is called the Eddington 
limit. 

The standard Eddington limit'^, was originally derived 
for stars. It imposes a fundamental upper limit for the 
radiative luminosity of a steady spherical accretion flow, 
which is usually called the Eddington luminosity. The 
basic idea is to find the equilibrium situation between 
infalling material when interacts with the outward radi- 
ation fiux. Assuming that the accreting material is fully 
ionized hydrogen, radiation pressure acts mainly on the 
free electrons. The associated force over one electron 
is due to Thomson scattering and it is given by aTF/c, 
where ar = 6.7 x 10^^'"' cm^ is the Thomson cross section, 
F the outward flux (in erg s~^ cm^) and c the speed of 
light. The radiation pressure pushed up electrons against 
the gravitational force, at the same time electrons drag 
the protons with them due to Coulomb interaction, but 
the force of radiation into protons is negligible next to 
electrons because the scattering-cross section of protons 
is [me/rnpf ~ 25 x 10~^ smaller than for electrons. 

Under this condition, accretion is only possible if the 
total gravitational attraction overcomes the radiation 
outward force. The Eddington limit corresponds to the 
situation when these quantities are equal. The gravita- 
tional force exercised from the central object of mass M 
to a electron-proton pair is given by GM{me + mp)/r'^ ~ 



^The standard or classical Eddington limit, introduced by Sir 
Arthur Stanley Eddington, only takes into account electron scat- 
tering processes. 



GMmp/r^ , and must be equal to the radiative force over 
one electron. Then we have 



GMtrip _ axF 



(42) 



The flux F of the accreting source is related to his lu- 
minosity through, F = L/Airr^ (as mentioned above, the 
system is spherically symmetric), then the last equation 
becomes 



GMm„ 
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c Airr^ 



^0. 



(43) 



The luminosity that fulfills this condition is the Ed- 
dington luminosity. 
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\.2Qxl{f^{M/MQ) ergs"^ (44) 



If we consider that an object (black hole, white dwarf, 
neutron star, etc.) is only powered by (spherical) accre- 
tion, this condition implies a limit on the steady accretion 
rate, called the Eddington rate Me- It can be obtained 
from the Eddington luminosity. 
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(45) 



When accretion exceeds this limit, the associated lu- 
minosity will exceed Le^ then the outward pressure radi- 
ation overtakes gravitational attraction hence; accretion 
is halted, turning off the source. 

For realistic astrophysical models, this approach needs 
to be reconsidered. Several of the used assumptions could 
not be true. For example, a mixture of elements other 
than hydrogen, or the ionization level of the gas invali- 
dates the latter. Spherical accretion seems more difficult 
to be produced. Viscous gas tends to form a disk when 
falling into a potential well due to loss of angular mo- 
mentum, therefore, the geometry ceases to be spherical. 

In the case of accretion disks, the Eddington limit 
should not be used as exposed in this section. 

2.3. Temperature profile and disk continuum 
spectra in a-disks 

Viscous processes in the internal region of the disk 
(equatorial plane) converts the potential energy of the 
accreting matter into thermal energy which is then com- 
pletely radiated away throughout the surface of the disk. 

Viscous heating per unit volume is given by: 



pv r 



dn 
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dr 



(46) 



where VIk — {GM/r^Y/'^ is the kcplcrian angular veloc- 
ity. 

The radiative flux F in the z-direction is given by the 
diffusion equation -which is only valid for optically thick 
disks-, therefore we have. 



oz 



(47) 



where D — Ac/3 is the diffusion coefficient, A the photon 
mean free path and u = aT^ correspond to the radiation 
energy density (a is the radiation constant). Notice that 
from the radiation energy density, du ~ AaT^dT, there- 
fore we can rewrite the above equation as (transforming 
the partial derivative into an ordinary one): 



F = -D 



= -D 



. du 

dz 
du dT 



dT dz 
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4XacT^ dT 
3 dz 
AacT^ dT 
Snp dz ' 



(48) 



where we have used for the photon mean free path 
A = 1/ up, with K. the Rosseland-mean opacity. For high 
temperatures (T > 10^ K) the opacity main sources are 
electron scattering and free- free (or Bremsstrahlung ^) 
absorption (for low metal abundance bound-free should 
be included). 



K= Ks + Kff = Ks + KopTj^^, 

0.4 cm^ g~i and kq = 0.64 x 10^^ 



(49) 



where «„ — U.4 cm' g " and kq — U.t)4 x iU'" m cgs 
units. The opacity of the disk is related to the optical 
depth through. 



kE 



(50) 



Now, from the energy balance equation, the radia- 
tion ffux must be compensated by the viscous dissipation 
heat, then we have 



■^German word, from bremsen "to brake" and strahlung "radiation" 
i.e. "braking radiation" or "deceleration radiation") 



dF 

dz 



(51) 



Integrating the above equation along the z-axis, we 
simply have 



rpi 



iQdis 



(58) 



Now from equations (54) and (56), we obtain the ef- 
fective temperature Teff{r), which is given by: 



F 



H/2 



q'^dz, 



(52) 



with the boundary conditions F{z = 0) = 0, and 
F{H/2) = aT^ff = q^H/2 (assuming black body emit- 
tance), where a is the Stephan-Boltzmann constant, Tg// 
the effective temperature, and H the total scale of height 
of the disk. 

From equations (48) and (52), we obtain the required 
equation for the temperature of the disk: 



AacT^ dT 
Sup dz 



= q' z. 



(53) 



The above equation means that a local energy balance 
between dissipated and radiated energy is assumed, i.e., 

^diss ^rad' 

The vertically integrated heat rate per unit surface 
Qdiss' ^^ obtained from. 



Q 



dis 



/ q+dz = -V^^], = -T^Tr^^K, (54) 



where the viscous tensor T^p (in a steady situation), is 
given by. 
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(55) 



Let us recall that in this prescription the term vT, is 
written as: 



The radiative cooling is given by; 



2 X F 
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(56) 



(57) 



where the factor 2 represents radiation from the two sur- 
face of the disk. Notice that the radiative flux of the 
disk scales as ex T^rr, a typical characteristic of black- 
body emissions. Using the above equation, along with 



the condition that Q^^^^ 
five temperature: 



ad' 



we obtain for the effec- 



T^fA^) 
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(59) 



For regions far away from the centre (i.e., r 3> riso), 
the effective temperature scales as: 



T, 



eff 



„-3/4 



(60) 



The radial dependencies of the effective temperature 
shows that the surface of the disk is hotter in the inner 
regions, and the disk is cooled as one moves away from 
the centre. Notice also that the temperature at the last 
stable orbit is equal to, Teff{riso) ~ 0. 

2.3.1. Thomson scattering case 

If we assume that the opacity is only due by Thomson 
scattering (k ~ 0.4 cni^ g^^), which is a good approx- 
imation for highly ionized disks (T ^ 2 x IC* K), the 
temperature profile is obtained by integrating Equation 
(53). Using the relation ac — 4cr and the boundary con- 
ditions mentioned before, after a simply calculation we 
obtain: 



T ^T, 



eff 



l + ?Wl-4±5 



-,1/4 



(61) 



The temperature at the equatorial plane is obtained 
taking z = in the above equation, then, 

/ 3X1/4 
Tc = T,ff (^1 + -r j , (62) 

and the optical depth of the disk is given in this case by. 



T = Ts = p— = S, 

niH 2 2TO/f 



(63) 



where cry is the Thomson cross section, and mn the pro- 
ton mass. 

Now, replacing the effective temperature Teff{r) (Eq. 
59) in Equation (61) or (62) we obtain the temperature 
profile or the temperature at the equatorial plane of the 
disk, respectively. 
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2.3.2. Disk spectrum 

We have already discussed that in the classic theory 
of accretion, the disk is assumed to be optically thick, 
radiating throughout the surface as a blackbody. The 
blackbody radiation is given by the Planck function. 



B.{T{r)) 



2h 



c^ eyiY){hv / ksTef f{r)} - 1 



(64) 



The flux is calculated locally, so at each different ra- 
dius a different value for the temperature is given. 

The spectrum of the disk is then calculated by inte- 
grating the specific intensity B,^{T{r)) over the surface 
of the disk, hence: 



black hole binary) shows soft and hard components (Mit- 
suda et al. 1984). The temperature associated to the soft 
part can be well fitted within this model, as shown below. 

Mineshige et al. (1994) fitted the soft X-ray spec- 
tra of Nova Muscae assuming blackbody emissions with 
a temperature profile T ex r~". Deriving the tempera- 
ture gradient from the X-ray data source, they show that 
— 91ogT/91ogr = n « 0.75 supporting the idea that soft 
components of the spectra are emitted by a classic accre- 
tion disk (see Fig. 2). 

The disk in Nova Muscae is not steady, but its vari- 
ability is slow enough so that time variations do not affect 
the emitted spectra, in fact Cheng et al. (1992) found 
that the UV-optical spectrum emissions are well fitted 
assuming a steady disk model. 



B,{T,ffir))2Trrdr, 



(65) 



where riso is the inner boundary of the disk (last stable 
orbit radius) and r^ his maximum extension. 

The effective temperature can be assumed to scale as 
a power-law function: 



Teffir)-To{^)\ 



(66) 



where To and rg are some normalization values. 

Replacing the above expression in Eq. (65), and using 
the Planck function B^{T{r)) (Eq. 64), we obtain for the 
spectral luminosity of the disk: 



L,. 



r{ 
n 



^Wo'/"/ B,{T)T-''^/''^-^dT, 



riso 



(67) 



wereXo = (^)(4)(f^)^/" 

From the above, we notice that the continuum spec- 
trum scales as: 



Li, oc V' 



,3-2/n 



(68) 



Using the theoretical value for the slope n of the tem- 
perature in the standard model {n — 3/4) (Eq. 60), we 
get for the spectrum L^, ex v'^/'^. Notice that the calcu- 
lated spectra is only valid in regions far away from the 
central object, where lower temperatures produce low en- 
ergy photons (i.e., soft X-rays, UV, optical). 

These results can be applied to several observations. 
For instance. X-rays spectra from Nova Muscae 1991 (a 






B 



1.0 
0.9 



0.6 - 

0.5 

O 



+ ^U 



30 1 00 

Total Uay in 1991 



0.90 
T^ 0.80 - 

bo 

O 






0.70 



'V O.60 



~T— r— r— r- 



+H-f^t 



50 1 00 

Total Oay in 1991 



Fig. 2. — Spectral fitting of Nova Muscae 1991 with the 
standard disk blackbody model. The top shows the disk 
temperature variability. The bottom the temperature 
exponent. Credits: Mineshige et al. (1994). 

The bolometric luminosity can be computed from the 
dissipated energy rate per unit volume q^ , therefore: 



L 



rd r^d 



q+dV = 27r / q+rHdr 



(69) 



Using Eq. (46) and (56), the above equation can be 
recast as: 
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L = 



3GMM 
2r2 



1 - 



ris 



GMM 



The above equation indicates that the disk luminosity 
is proportional to the accretion rate of the black hole; 
higher rates, produce more luminous disks. Notice also 
that the factor GM /2riso correspond to the half of the 
potential energy difference between infinity to the last 
stable orbit. 

Since the accretion rate originates from turbulent vis- 
cous processes, the ability to convert potential energy 
into radiation, comes exclusively from viscous work. 

It is worthy to note that, since the disk is hotter in 
the inner regions, and the spectra of the disk show es- 
sentially blackbody emissions, we conclude that in this 
model, high energy photons are emitted from the inner 
regions whereas low energy photons come from the outer 
portions of the disk. 

2.4. Uncertainties in a-disk models 

As we have seen, Shakura-Sunyaev formalism explains 
reasonably well disks present in binary systems but such 
a model is faced to several stability problems and some 
assumptions must be revisited. 

We begin the analysis by identifying the typical timescales 
on which the disk structure may vary. 

The shortest one corresponds to the dynamical timescale, 
which corresponds to the orbital timescale 



to == ^' 



(71) 



where fi corresponds to the angular velocity, which is de- 
duced from the gravitational potential ^ from il^ — \^- 
For a Keplerian, non self-gravitating disk, the angular 
velocity is given by, il^ = Q.\ = ^^, where M is the 
mass of the central object. In this case, the dynamical 
timescale is simply given by. 



to 



m' 



(72) 



which coincide with the timescale needed for inhonio- 
geneities to be dissipated (be smoothed) in the z-dircction. 
In other words, deviations from hydrostatic equilibrium 
in the vertical direction are smoothed in a timescale nec- 
essary to a sound wave with speed Cg , to scroll a distance 
H. This time is given by. 



_ H _ H _„_i _ 

tr — — TT — i* — tjl-, 



(73) 



where we have used the fact that in a thin standard disk, 
the sound speed is related to the scale of height as, Cg — 
Hfl. 

The dynamical timescale can be compared with the 
viscous timescale in order to estimate the effect of vis- 
cosity on the evolution of accretion disk. The viscous 
timescale gives the timescales in which matter diffuses 
through the disk over a distance r, under the effect of 
viscous torques. This quantity is given by. 



r 

V 



inserting the a-ansatz in this equation, we obtain 



til 



ac^H 



(74) 



(75) 



Another important quantity is the thermal timescale, 
which gives the characteristic time needed for the system 
to recover thermal equilibrium (when Q^ = Q~). The 
thermal timescale is given by a ratio between the thermal 
energy content (^ pkT/fxmH ^ pc^) per unit surface area 
and the viscous dissipation rate {Q~^ « rHTr,pdQ/dr « 
{9 / 4)il'^ h'Y,) per unit surface area, then. 



Hh 



4- = M-\ 



(76) 



where M — v^pj Cs is the Mach number and t^ the viscous 
timescale given by Equation (74). From the condition im- 
posed in this model, that the viscous energy dissipation 
is locally balanced by radiation from the two disk sur- 
faces, the cooling timescale, which is given by the ratio 
of the thermal energy content by the radiative lost rate 
(Q" — aT^ / KfipH'^, when the disk is optically thick) 
should be the same as the thermal timescale. When ad- 
vection dominates the disk, this equality does not hold, 
and the cooling time will be greater than the thermal 
timescale. 

Now, we can obtain some relations between these 
characteristic timescales. We have already found that 
to = tz, and tth ~ M^'^ty. From the viscous timescale 
(Eq. 75), we can rewrite t^ « —^—tf — ■r^-r-'^, then, 



aCsH 



Summarizing the results, we see that, 
Ito^tz^ atth « a f — j ty 



(77) 



(78) 



and considering that a ^ 1, the hierarchy for numerical 
values is 
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to ~ tz <, Uh ^ it)- 



(79) 



It is common to see astrophysical m.odels based on 
steady accretion a-disks. Generally, these models teach 
us important features of these objects. A steady solu- 
tion means that the system is, in fact, in an equilibrium 
situation, and must stay unchanged under small pertur- 
bations; otherwise, the steady solution cannot exist in 
reality. 

We will see briefly two cases of instability situation 
in models based on Shakura-Sunyaev formalism, namely 
the thermal and viscous instability. 

The thermal instability occurs when the energy bal- 
ance between radiative cooling and viscous dissipation is 
perturbed, and this perturbation grows on a timescale tth 
much less than the viscous timescale t^, they take place 
specially in the inner region of (standard) accretion disk, 
where temperatures are higher. More clearly, suppose 
that the disk is in local energy balance (Q"*" = Q~) and 
we increase the central temperature Tc by a small pertur- 
bation 5Tc, what happens in this model is that, after the 
perturbation, the local heating and cooling rates increase 
(which is normal), but the local heating rate increases 
much faster than the cooling rate. This will lead to a 
thermal runaway, making the steady situation unrealis- 
tic. 

To easily see how these heating and cooling rates 
change differently, we consider a simple example: Sup- 
pose an optically thick and geometrically thin accre- 
tion disk (i.e., the a— model), where the cooling rate 
is given by Q~ oc T^ and the dissipation rate is given by 
Q+ ^ rHTr^dfl/dr. Considering the parameterization 
Tr^ — —aP (Eq. 40) and a radiation pressure P ex T**, 
in a vertical hydrostatic equilibrium with H ^ 2P/f2if S, 
we obtain for the heating rate Q+ ex T*, while for the 
cooling rate Q~ ex T^. Comparing these two expres- 
sions, we notice that when T increases, the cooling rate 
will grow more slowly than Q^ . Hence, any perturbation 
in the temperature will lead to a growing instability. 

Unless very specific and often unrealistic conditions 
are supposed, instabilities cannot be avoided. 

The general condition for thermal stability can be 
written as. 
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(80) 
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whenever the above inequality is satisfied the system 
shows a thermal instability. 

If we consider changes in the disk structure, it can be 
shown that a perturbation ST, in the steady solution Eq 



(i.e., S = Eq + <5S)i the the viscosity perturbation 6fi 
(obtained from /i = i^E), obeys a diffusion equation, and 
whenever the condition dii/dT > is locally satisfied, 
mass will be removed from that region, leaving a frag- 
mentation of the disk breaking up in a series of rings. A 
general condition for disk fragmentation is given by. 
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(81) 



where M is the accretion rate. In standard accretion 
models, the relation M = M(S) takes the S-shape, which 
indicates an instability region whenever the slope of the 
curve is negative (see Figure 3). 
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Fig. 3. — Schematic S-shape plot with the instability 
occurring on the intermediate branch of negative slope. 

A similar shape emerges in the analysis of T = T(E), 
showing the development of thermal instabilities starting 
from relatively low temperatures (e.g. T > 10^ K). 

For a complete stability analysis on Shakura-Sunyaev 
based models, see Piran (1978), and Pringle (1981). 

Conclusions arising from stability analysis, is that 
steady-accretion disk with the a— prescription cannot be 
really stationary. Another uncertainty in these models 
is about the type of pressure that should be used in the 
anomalous stress tensor (Eq. 40); It should be used only 
the gas pressure? radiation pressure? something else? 
Even just considering gas pressure, and hence avoid- 
ing thermal and viscous instabilities, radiation pressure 
is still important. In fact, in Shakura-Sunyaev models 
around black holes, radiation pressure is always domi- 
nant for high accretion rates (when the luminosity is near 
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the Eddington luminosity), so it should be included, but 
their inclusion leads to thermal instabilities. 

2.5. Beyond the standard model 

We have already seen some drawbacks in the a- 
parameterization. It is clear now that a key ingredient 
to describe the accretion process is the turbulent viscos- 
ity. It causes the angular momentum loss and the inward 
motion of the material (as well as an outward flow by the 
external region of the disk). 

Shakura-Sunyaev parameterization remains an ad-hoc 
description limited to thin disk geometry and negligible 
disk masses. Now, if you want to "construct" a supermas- 
sive black hole through an accretion disk in timescales 
less than 1 Gyr, needed for consistence with observa- 
tions of most distant quasars at redshifts z ^ 6, you will 
need to consider a self-gravitating initial disk. We ex- 
pect from this a high accretion rate, that could lead to 
a super-Eddington regime, and probably an inner region 
with disk thickness H{r) comparable to the disk size i.e., 
H ^r. 

For the reasons mentioned above, we need a more real- 
istic description of the turbulent nature of viscosity, and 
relax the assumption that the disk can only be geomet- 
rically thin. 

Understanding the origin of turbulent transport is one 
of the major problem in accretion physics. The nature of 
the physical mechanism which causes turbulent motion 
is still an open question. A good candidate that partially 
explain turbulent stress is the magneto-hydrodynamic in- 
stability (MHD instability). It was originally discovered 
by Velikhov (1959). From this, several authors explored 
this possibility, and applied they results to accretion disk 
theory. 

In these models, the source of turbulent viscosity may 
arise from magnetic fields, which are generated and main- 
tained in the disk by a dynamo mechanism. This mech- 
anism could produce local instabilities generating tur- 
bulence and hence an a— effect (with a magnetic pres- 
sure), similar to those produced by the standard model 
(see for instance Brandenburg (1995), Balbus & Hawley 
(1998). Numerical simulations showed that a weak mag- 
netic field in accretion disks could lead to a very effective 
transport of angular momentum and hence, an efficient 
accretion rate for a sufficiently ionized disk fiow (e.g., 
Hawley (1995)). 

These models work quite well, they describe very eas- 
ily low-mass X-ray binaries, or systems when the cen- 
tral object is actually a magnetic neutron star or a non- 
magnetic simple star. It is important to note that in 
those cases the disk is probably heated by irradiation 



due the central object, ensuring a sufficiently ionization 
rate. 

The situation is rather different when a black hole is 
considered as the central source, and even more radical 
if we attempt to describe large scale disks (^ hundreds 
parsec radius), where it is expected that the outer region 
of the disk the temperature is sufficiently low to form a 
large amount of molecular gas. Then, we cannot always 
expect a high ionization rate in the gas and therefore the 
necessary conditions to anchor a magnetic field in it. 

On the other hand, due to the turbulent nature of 
the flow, it is expected that the number of Reynolds be 
relatively high. This number is deflned as. 
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(82) 



where u is the characteristic mean fluid velocity, / the 
characteristic dimension, and v the kinematical viscos- 
ity (related to the dynamical viscosity through v — ^/ p 
measured in m^ s^^). 

Lyndell-Bell and Pringle (1974) pointed out that if 
an accretion disk is driven by molecular viscosity (where 
v — )■ small) , the Reynolds number becomes high (see Eq. 
82), then turbulent eddies develop in the flow, enhanc- 
ing in turn the effective viscosity (turbulent viscosity), as 
a consequence the Reynolds number will decrease. This 
process will continue until a critical Reynolds number iRc 
is reached, which indicates the onset of turbulence. In 
this sense, the flow is expected to be self-regulated and 
characterized by a critical Reynolds number. Below 3ffc 
the flow is laminar, but for higher values the flow becomes 
turbulent. This ensures that if the angular momentum 
transfer is inefficient, turbulence decays triggering insta- 
bilities which regenerates the turbulence. 

This approach was revisited by de Freitas Pacheco & 
Steiner (1976 ), who modelled a steady accretion disk 
able to explain the observed X-ray spectrum of the binary 
source Cygnus X-1. Their prescription is robust, assuring 
high accretion rates and making magnetic fields needles. 

Notice that a similar treatment was discussed by 
Duschl, Strittmatter & Biermann (1998) (see also, Richard 
& Zahn (1999)), who have introduced the "/3- viscosity" , 
which has essentially the same parametrization as that 
of de Freitas Pacheco & Steiner (1976). 

The best connexion with experimental data comes 
from laboratory experiments on rotating Couette- Taylor 
flow (Wendt 1933), Taylor 1936), consistent in a fluid be- 
tween two cylinders rotating at different speeds aiming 
to study hydrodynamical instabilities in the flow. 

Considering the Reynolds number definition (Eq. 82), 
we assume that the characteristic mean fluid velocity is 
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very close to the azimuthal velocity u "^ v^p and the char- 
acteristic dimension of the system is I ~ 27rr. After the 
system has adjusted itself to a critical Reynolds number 
5Rc, we get, 



27r 



•■^Lfl 



/3nr^, 



(83) 



where H. = v^/r corresponds to the angular velocity, and 
the free parameter (3 = 2tt/^c gives a measure of how 
turbulent the fluid will be. 

Based on hydrodynamical experiments, common val- 
ues for the critical Reynolds number are 5ftc ^ 10^ — 10*. 
Recall that the smaller 5Rc is, the earlier the turbulence is 
generated, and then more turbulent the fluid is, ensuring 
efficient angular momentum transport outside the disk, 
and allowing material to flow in. 

This parametrization of viscosity differs in many as- 
pects from Shakura-Sunyaev models. It allows an effi- 
cient transport of angular momentum like his predeces- 
sor, but under the analysis of stability given by Piran 
(1978), it shows a stable disk for all /3 values and for sev- 
eral cooling mechanism. For instance, see Montesinos & 
de Freitas Pacheco (2011), were using this parametriza- 
tion they studied the evolution of self gravitating disks 
around supermassive black holes solving numerically the 
complete set of hydrodynamical equations discussed in 
this paper. 

3. Evidence of self-gravitating accretion disks 

In the Universe there are several classes of self-gravitating 
disks. They can be small circumplanetary disks with less 
than 1 AU radius, or protostellar disks with intermedi- 
ate radius of (100 — 1000 AU), or several hundred parsec 
scale disk as found in AGNs. In all of them, self-gravity 
plays an important role in the physical properties of these 
object. 

After the study of 47 objects, Hillenbrand et al. (1992) 
argued that Herbig Ae/Be stars* are surrounded by cir- 
cumstellar accretion disks with masses up to six times 
higher than the mass of the central object. 

As another example of self-gravitating circumstellar 
accretion disk, Chini et al. (2004) reports observations of 
a massive protostar of about ^ 2OM0 fed by an accretion 
disk of about ~ lOOM© (see Fig. 4). 

We are interested in accretion disks in AGNs. Obser- 
vational evidence of these objects came from H2O maser 
emission at radius of about ~ 1 pc from the central black 
hole. From these emissions, the rotational velocity curve 




^A young star {< 10® yr) of spectral type A or B, with mass between 
1.5 to 10 Mq 
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Fig. 4. — Silhouette of the accretion disk with a diameter 
of 20 000 AU, seen against the light from the H II region. 
White contours delineate the densest part of the disk 
(inner torus). Credits: Chini et al. (2004). 



of the disk can be deduced, which differs from Keplerian 
rotation. 

A good example is the case of NGC 3079, classified as 
a Seyfert 2 (Ford et al. (1986), Ho et al. (1997) (Figure 5 
shows a picture of the whole galaxy) . From water maser 
emissions, the obtained disk rotation curve is relatively 
flat (Sofue 1997). They also indicates a central compact 
object with a mass of about ~ 2 x lO^Af© enclosed within 
a radius of ~ 0.4 pc. CO(l-O) observations suggest that 
the disk surrounding the black hole is relatively thick; 
with an aspect ratio in the range of H/r ^ 0.1 — 0.5, and 
a mass of the parsec-scale disk of about ^ 3 x 1O^M0, 
therefore self-gravitating (Koda et al 2002). 

Based on stability, cooling, and timescale considera- 
tions, Kondratko et al. (2005) also argue that the disk 
is self-gravitating and clumpy with the necessary condi- 
tions for star formation. 

A similar situation is observed in NGC 1068, another 
active galaxy classified as a Seyfert 2 that harbor a 
SMBH. From water maser emission observations. Green- 
hill & Gwinn (1997) show a non Keplerian rotation disk 
curve. The velocity field indicates that the rotational 
velocity appears to rise steeply in the center toward the 
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Fig. 5.— Galaxy NGC 3079. The structure is more than 
3000 hght- years wide and rises 3500 hght-years above the 
galaxy's disk, he smaller photo on the right is a close-up 
view of the bubble. Astronomers suspect that the bubble 
is being blown by winds released during a burst of star 
formation. Source: Hubble Space Telescope. 



peak on the molecular ring (Fig. 6), Lodato & Bertin 
(2003) provide an interpretation of the non-Keplerian 
rotation in NGC 1068 by means of a self-gravitating ac- 
cretion disk model where the resulting fit-curve is not a 
power-law (Bertin & Lodato 1999). They also derive a 
black hole mass of (8.0 ± 0.3) x IO^Mq. 

More recent observations from the VLTI suggest that 
the central SMBH is surrounded by a torus of gas and 
dust of size ^ 3.4 pc (Jaffe et al. 2004). It is expected 
that this torus-like structure should be self-gravitating. 

Self-gravity not only affects the rotational curve of 
the disk, it also influences its emitted spectrum. For 
instance, T-Tauri stars'^ are objects believed to be sur- 
rounded by a disk of gas and dust. The evidence for such 
disks comes from the excess of strong infrared (Rydgren, 
Strom & Strom 1976) and millimetric radiations (Beck- 
with et al. 1990), indicating that these emissions come 
from regions far from the central star. 

The light curve in T-Tauri stars can be described by 
a power law model in the form L^ ex ly'^'^ (Adams et 
al. 1988). If we take a disk, whose temperature follows 
a power law T ex r~", the associated luminosity takes 
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®The last stage of stellar evolution before the main sequence, they 
are characterized by emission lines, rapid variability, and X-ray 
emission, with masses of roughly 0.2 to 5 Mq 
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Fig. 6. — Fit to the rotation curve from water maser 
emission by a self-gravitating accretion disk model. The 
small panel shows a blow up of the declining part of the 
rotation curve together with the best fit obtained by as- 
suming Keplerian rotation. Credits: Lodato & Bertin 
(2003). 



the form Ly ex v'^~'^l°' (see Section 2.3), therefore the 
temperature must decrease with the radius in the form 
r~^/^ to be in agreement with observations of the spec- 
trum. The problem arises when assuming a standard 
a-disk, which implies a Keplerian rotation curve. The 
temperature profile is obtained when the viscous dissi- 
pation equals the radiative dissipation i.e., Q\igg = Q^ad 
(assuming black body emissions), making the tempera- 
ture decrease as r~'^'^ contrary to observations. In or- 
der to explain T-Tauri spectrum, we need another "non- 
standard" disk model. A deviation from standard models 
that fit observational data quite well, is to assume that 
the rotational curve of the disk is not Keplerian. For 
instance, a T-Tauri star with flat spectrum, could be ex- 
plained with a constant rotational curve. This kind of 
shape appears naturally in Mestel's self-gravitating disk 
models (Mestel 1963). Following this reasoning, the in- 
frared excess could be a consequence of the self-gravity 
in the disk. Assuming a marginally stable disk, with 
Toomre parameter Q — 1, gravitational instabilities pro- 
vide a supplementary contribution to the heating of the 
disk, which coincides with the infrared excess (see for in- 
stance Bertin & Lodato (1999), Lodato & Bertin (2001)). 
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The spectral properties of AGNs are niodified in the 
same way when self-gravity is taken into account. Kawaguchi 
et al. (2004), explain the origin of Mid- Infrared to X-ray 
spectrum emissions of the Narrow-Line Seyfert I galaxy 
TONS 180, assuming that the disk can be divided in 
three different zones. An inner region of about ~ SOOOr^ 
from the black hole, where the disk is slim^. An interme- 
diate self-gravitating region, where gravitational instabil- 
ities and the non-Keplerian rotational velocities play an 
important contribution to viscous heat, and an external 
region composed of a dusty torus-like structure, probably 
at the origin of infrared radiation. Their considerations 
show a very good fit with observations (see Figure 7) , es- 
tablishing clearly the importance of self-gravity in AGNs 
models. 



45 



43 



Log(hv[keV]) 
-3 -a -1 



-¥-¥- 



IR 



•^ 



-^ I ^ ^ < i I I 



1- 
Optical 



-a 







\ \ \ -i \ 



X-ray 



-+i-+ 



-SG disc 
corona 



SG disc 



■^ Standard 'ife 



_i I I i_ 



lO/xm ■ Ifim lOOOA 



disc 



H 



''^i 



% 



13 



14 



15 16 

Log(i/[Hz]) 



17 



18 



Fig. 7. — Mid-IR to X-ray spectral modeling of Ton S 
180. The open (white) squares and thick striped lines 
are the observed data. The two solid lines (green and 
magenta) are the spectral components of the dusty torus 
(left) and the non self-gravitating disk. The total spectra 
of those three regions are indicated by dashed lines: the 
thick line (red) is the spectrum with the self-gravitating 
disk, while the case with a standard disc is indicated by 
a lower curve. Credits: Kawaguchi et al. (2004). 



3.1. Self-gravitating disk model 

As we have seen, self-gravity plays an important role 
in the evolution of astrophysical disks, however not many 
simulations in accretion disk include it. This could be due 
to the difficulty to solve properly the Poisson equation or 
the belief that self-gravity could not be important if the 



^Accretion disk model introduced by Abramowicz et al. (1988) 



mass of the disk Aid is of the same order or just a few 
percent higher than the mass of the central object M,. 
This last can be true if the disk is very light, but actually 
it is not the ratio between the mass of the central object 
and the disk mass that needs to be compared in order to 
include or exclude the self-gravity from a model, but the 
ratio of their gravitational potential. 

At the present time, there is no simple analytical ex- 
pression for the gravitational field created by a disk. The 
best tools are numerical calculations. 

We need a mathematical model for the self-gravity of 
a disk. Therefore, we need to solve Poisson's equation. 
They are several ways to solve this equation in order to 
compute the gravitational field ^(x) created by a mass 
distribution of density p(x), where x is a space- vector. 
But this does not mean that the task is easy. It is spe- 
cially hard in astrophysical problems, where the matter 
distribution is not completely known. This is why, sev- 
eral assumptions, based on observations, must be consid- 
ered. 

Generally speaking, Poisson's equations is, 

V^*(x) =47rGp(x), (84) 

where ^(x) its the gravitational potential created by 
/9(x), and G is the gravitational constant. 

The most general solution for this equation takes the 
form, 



*(x) 



G 



p(x')d3x' 



V 



(85) 



and the integration is computed within the volume V 
enclosing the density distribution p(x'). 

The acceleration g(x) produced by the gravitational 
potential 4'(x), can be computed from the relation g(x) = 
— V*(x). Then we have. 



g(x) = -G f 



p(x')(x-x')d3x' 



(86) 



Using symmetry arguments, and galaxy observations 
(like observed rotational curves) we can assume, for sim- 
plicity, an axisymmetric model (i.e., 3/ dip = 0). This im- 
plies that the above space- vector x, refers to x = (r, (p^ z). 

Next, we consider the gravitational potential of an ho- 
mogeneous spheroid with columnar density profile S(r) 
and then make the necessary limits and extrapolations to 
obtain the potential on a flat disk approximation. These 
calculations are based on the work of Wyse and Mayall 
(1942); E.M. Burbidge, G.R. Burbidge and Prendergast 
(1959), whose purpose was to study the relation between 
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the density distribution on disks (or geometries close to 
it) and their gravitational field, aiming to reconstruct ro- 
tational curves in the nuclear region (inner region where 
dark matter plays no appreciable role r < 3 Kpc, see Fig- 
ure 8) from several observed galaxies like, M31, M33, spi- 
ral nebula NGC 2146, NGC 3198, amount others. Their 
analysis gives a recipe on how to calculate Newtonian 
gravitational potentials produced by some density distri- 
butions in a disk geometry. 



DISTRIBUTION OF DARK MATTER IN NGC 3198 
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Fig. 8.— The rotation curve for the galaxy NGC3198. 
Note that in the inner region, observations fits very well 
with the disk model. 

We start our analysis assuming an oblate spheroid 
(Fig. 9) with major and minor radii a and b respec- 
tively, eccentricity e = (1 — &^/a^)^/^ and a density p{r). 

In a radial equilibrium situation, centrifugal forces 
given by — u^/r, equal the gravitational pull F(r). Ne- 
glecting viscous forces and gradient pressure, we simply 
have. 




Fig. 9. — Oblate spheroid. 



^=F(r). 



(87) 



The total gravitational force F{r) at a point r can be 
expressed as. 



F{r) 



/'"<•■''§*■■ 



where the function ^ = C(r; a, e) can be obtained from 
the expression of the gravitational potential for an ex- 
terior point. In that region (where p{r) = 0), Poisson's 
equation reads d^'^/dr^ — 0. Therefore, we obtain: 



2 2 2n1/2 

r — e a ) ' 



■ -1 /eay 



/■I \ -^'"'^/i 2n1/2 

(89) 

From the above equation, we take the derivative dC,/dr, 
and we replace it in the equation for F{r) (Eq. 88). Us- 
ing the equilibrium Equation (87), we find that 



vl{r) 



47rG(l - e^) 



2N1/2 



p{a)a^da 

(^2 _ £2^2)1/2 ■ 



(90) 



In Eq. (90), the integration variable a, describes an 
equipotential satisfying the relation a^ = x^ + y'^, while 
r^ = x'^ + y'^ + z"^ . Taking now the limit for a flat disk 
approximation, which means that the ratio between the 
minor and major radius tends to zero, i.e., b/a -^ 0, and 
the volume density p{r) must be replaced by the surface 
density S(r). We obtain. 



vlir)=G 



Y,{a)a'^da 



(r2 



2)1/2 



(91) 



Using the relation d'^/dr = —v'i/r we express the 
gravitational potential created by the disk, as; 



dr 



G 

r 



J:ia)a^ 



(r2-a2)i/2 



dr. 



(92) 



This is the well know flat disk approximation, for fur- 
ther reading about this approximation see Brandt & Bel- 
ton (1962), and for mathematical details about the po- 
tential theory O.D. Kellogg (1956). 

From a direct calculation, we can obtain the mass dis- 
tribution of the disk M^ir) inside the radius r by inte- 
gration of the differential mass dM{r) — 2TTr'S{r)dr, thus 
obtaining 



Mdir) 



2TT'E{a)ada. 



(93) 
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We calculate the radial acceleration due the gravita- 
tional influence of the disk (i.e., g^ = —d'^ disk/dr) using 
the Eq. 92. As for the vertical component of the gravi- 
tational acceleration produced by the disk, we have 



disk 



9z 



= 2^GE-, 



(94) 



where H is the effective scale of height of the disk. 

The contribution from the black hole to the radial and 
vertical acceleration, g'^'^ and gl'^ are, respectively. 



^hh 



9r 



GMkh 

(r- rs)2 



and 



^bh 



GMbh 



(95) 

where rg — 2GMi,h/c^ is the Schwarzschild radius. The 
expression for 5^'' comes from the Paczynski-Wiita Po- 
tential Paczyhski & Wiita (1980). 

From these expressions, a simple parameter can be 
used to quantify the influence of the self-gravity. We 
define two quantities by the ratios; ^r = ffr^^'^/ffr'') ^"^^ 

^ — ^disk / ^,bh 
'^z — J/z /i/z • 

If we consider an initial surface density profile given 
by a power law, 



S(r) = Eo 



(96) 



the expression for (^z can be greatly simplified: 

Firstly, we note that from the equation for the mass 
of the disk (93), we have 



Ma{r) - / 27rEo 



dr = 27rEo 



ro 



2-n 
(97) 



The expression for <;z now reads, 



27rE(r)r3 



27rEo 



HMbh 
and using Eq. (97), we obtain 



\HJ Ml 



Mbh' 



(2- 






(98) 



(99) 



The above equation indicate that the disk becomes 
self-gravitating in the vertical direction when Md{r)/Mbh {r) 
H/r. Under normal conditions, the aspect ratio H/r of 
the disk takes values in the range ^ 10^^ — 10^^, indicat- 
ing that even when the mass of the disk is smaller than 
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the black hole mass, self gravity plays an important role 
and should be considered. 

A similar condition is obtained in terms of the Toomre 
parameter Q (Toomre 1964). Toomre established the 
basic principles concerning the gravitational instabilities 
in accretion disks. A gaseous disk is locally unstable to 
axisymmetric perturbations, if 



Q 



kCg 

ttGD 



<1, 



(100) 



where k is the epicycle frequency. This frequency comes 
from small oscillations of a particle motion (infinitesimal 
perturbations) in a stable circular orbit without chang- 
ing its angular momentum. The radial motion around the 
central object is periodic in that orbit, called epicyclic os- 
cillation, his frequency is then called epicyclic frequency. 
It is defined by: 



k^ = 20M 2 



dlnn 
dlnr 



(101) 



where the angular velocity fl should be taken from the 
total gravitational potential ^t = 'i'bh + "^disk, using 
n^ — ^^S^. In the case of Keplerian rotation, we have 

fc - Ok = ./GMbh/r^. 

To see more clearly what lies behind the Q parameter, 
we take the following approximation: We assume that the 
disk approaches a Keplerian rotational velocity (which, of 
course, is not realistic for a self-gravitating regimen), i.e., 
17 PS CIk. In this case, the epicycle frequency is /c ~ VIk 
and Cs ~ Hflii . Now, from the definition of the Toomre 
parameter Q, using the density profile given by (96), and 
the mass of the disk Md{r) as a function of S(r) given in 
Eq. (97) (assuming n = 1), we have 



Q - 



ttGH^ 



r 
r 



ra 

Mbhjr) 
Md{r) 

Mbhjr) 
Md{r) 



(2-n) 



(102) 



Again appears the ratio Mbh{r)/Md{r), indicating 
that the self-gravity begins to be important starting from 
Md{r) / Mbh{r) ~ H/r. Recall that the above approxima- 
tion was setting assuming Keplerian motion, therefore 
k = Q.K] so, what Eq. (102) indicate is a limit between 
a self-gravitating and a Keplerian regime for the scale of 
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height. The Toonire parameter Q, defines the gravita- 
tional instability threshold under the condition Q ^ 1 for 
stable regime. 

3.2. Other methods 

There are also other methods to solve Poisson's equa- 
tion. For instance, development over the Green function, 
in which the term 1/(| x — x' |) in the Eq. (85) is decom- 
posed as. 



4-7r T 



(I x-x' I) ^-^ ^-^ 2/4-lrL 

(103) 
where r< = Inf (r,r'), r> = Sup(r, r') and Yim{0,(p) are 
the spherical harmonics used as base, see for instance 
Cohl & Tohline (1999). 

A similar approach can be done using Bessel functions 
or using Fast Fourier Transform (FFT) methods (e.g. 
Binney & Tremaine (1988), this last one is very pop- 
ular in 2D self-gravitating accretion disk models. The 
drawbacks in the FFT approach are linked to the preci- 
sion (is a first order method), but this can be improved 
increasing grid number. This is not an issue for relatively 
small gaseous disk models (as protoplanetary disks) , but 
for larger disks (ranging from km to pc) numerical cost 
using those methods are unavailable for most of the com- 
puter clusters. 
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